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High energy scattering in the QCD parton model was recently shown to be a reaction-diffusion 
process and, thus, to lie in the universality class of the stochastic Fisher-Kolmogorov-Petrovsky- 
Piscounov equation. We recall that the latter appears naturally in the context of the parton model. 
We provide a thorough numerical analysis of the mean field approximation, given in QCD by the 
Balitsky-Kovchegov equation. In the framework of a simple stochastic toy model that captures 
the relevant features of QCD, we discuss and illustrate the universal properties of such stochastic 
models. We investigate, in particular, the validity of the mean field approximation and how it is 
broken by fluctuations. We find that the mean field approximation is a good approximation in the 
initial stages of the evolution in rapidity. 

I. INTRODUCTION 

The description of the rise of hadronic cross sections with the center of mass energy ^fs within quantum chro- 
modynamics (QCD) has been a challenging issue. Old experimental results, e.g. on proton scattering, do not yet 
have a satisfactory theoretical explanation. So far, no realization has been found for the Froissart bound, which is 
a consequence of the unitarity of the .S'-matrix and which is a bound on the high energy behavior of the total cross 
sections, a < (1/to 2 ) In 2 s. 

A perturbative calculation of the evolution of scattering amplitudes with energy, achieved by Balitsky, Fadin, Kuraev 
and Lipatov (BFKL) 30 years ago [| and subsequently improved to next-to-leading order fj, is apparently successful 
in describing the rise of photon-proton cross sections when the virtuality of the photon is high enough. However, the 
BFKL equation, being linear, is not compatible with the Froissart bound when it is extrapolated to very high energies. 
Complying with the limits set by unitarity requires the introduction of nonlinear terms in the evolution equations. A 
first stepin this direction was taken by Gribov, Levin and Ryskin (GLR) in 1981 Q, and by Mueller and Qiu (MQ) 
in 1986 Q. Subsequently, more involved QCD evolution equations were derived within different formalisms. Balitsky, 
Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov and Kovner (B-JIMWLK) HH] have developed a comprehensive 
approach to high energy scattering from different but converging points of view. Technically, they were able to write 
the high energy evolution of QCD amplitudes in the form of a functional integro-differential equation, that can also 
be expressed as a Langevin equation or as an infinite hierarchy of coupled differential equations. A much simpler 
equation was obtained by Balitsky [j| and by Kovchegov (BK) in the particular physical context where the target 
is a large nucleus. The latter turns out to be a specific limit of the former, and its structure is similar to the GLR-MQ 
equation. 

Meanwhile, on the phenomenological side, Golec-Biernat and Wusthoff showed that unitarization effects, that 
cannot be taken into account by a linear evolution such as the BFKL equation, may have already be seen at the DESY 
electron-proton collider HERA. Their observation gave a new impetus to the field and triggered many phenomeno- 
logical and theoretical studies. Furthermore, the model they had proposed pointed towards a new scaling, "geometric 
scaling" , which was subsequently found in the HERA data [9j . 

Much insight and intuition has been gained from numerical studies: even before the B-JIMWLK formalism had 
been developed, Salam performed numerical calculations of scattering amplitudes near the unitarity limits in the 
framework of a Monte-Carlo implementation E3] °f Mueller's color dipole model E3- More recently, much effort 
has been devoted to obtaining numcr ical solutions of the BK and B-JIMWLK equations El El El El El, 

and in 

particular to the understanding of geometric scaling. But no analytical solution could be found. 

Important progress was accomplished when the first terms of the large rapidity asymptotic expansion of the solutions 
to the Balitsky-Kovchegov equation were computed by Mueller and Triantafyllopoulos 17]. This expansion was 
subsequently systematized with the discovery 0, IT^I that the BK equation lies in the universality class of the Fisher- 
Kolmogorov-Petrovsky-Piscounov (FKPP) equation |2(J,|2l|], and that geometric scaling was equivalent to the property 
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that the latter admits traveling wave solutions. The FKPP equation describes a vast class of physical phenomena, 
actively studied by mathematicians and statistical physicists. 

Meanwhile, Mueller and Shoshi made a first step in going beyond the BK approximation to high energy 
scattering in QCD. Through a sophisticated analysis of the unitarity and symmetry constraints on the QCD evolution, 
they were able to propose a correction to the saturation scale, induced by effects neglected in the BK evolution, and 
realized that geometric scaling was asymptotically broken. 

The complete equivalence between high energy QCD and statistical physics models of reaction-diffusion type, first 
understood by one of us (S.M.) |2jJ,|2J|, provided a new picture in which the nature of the effects neglected by the 
BK equation became transparent. The universality class of high energy QCD was identified as that of the stochastic 
Fisher-Kolmogorov-Petrovsky-Piscounov (sFKPP) equation [2||. The latter applies to a very wide class of physical 
phenomena, ranging from directed percolation to population growth. It has many applications, e.g. to biology, 
genetics, anthropology, fluid mechanics. Roughly speaking, most of the physical situations in which objects evolve by 
multiplying and diffusing, up to a certain limiting threshold, may be captured by an equation in the universality class 
of the sFKPP equation. 

The crucial point is that the QCD parton model has such dynamics, when viewed in a particular way that we 
will explain in this paper. Once this observation is made, all mathematical results can be transposed to the physical 
context of high energy QCD |24|. In particular, Mueller and Shoshi's results on the saturation scale are recovered. 
The breaking of geometric scaling, that they had understood qualitatively in their approach, is shown to be related to 
the statistical dispersion in the occupation numbers of partons in the hadronic Fock state. This new point of view on 
high energy QCD not only confirmed and brought new exact results on QCD scattering amplitudes: it also provided 
a physical understanding of the very nature and the role of fluctuations contained in the QCD evolution equations. 
Moreover, it helped to understand how to set up a systematics to go beyond the BK equation. It was recognized that 
the B- JIMWLK formalism was not complete [2(| H?) , thus confirming what was expected 0, 13 HH • The form of the 
new terms can be written down taking advantage of the correspondence to statistical physics. 

The goal of this paper is to provide numerical solutions of the evolution equations of QCD, in order to study 
accurately their universal properties. The outline goes as follows: in Sec. ^ we explain how the parton model 
is related to a reaction-diffusion process. ^From physical but rigorous arguments, we derive in a simple way the 
evolution of the parton densities with energy, reviewing and expanding upon Refs. |23l |24|. In particular, we show 
how the sFKPP equation appears. In Sec. 11111 we study numerically in detail a useful mean field approximation to 
the full evolution equations — that is, the BK equation — in relation to the most recent analytical results. In the final 
Sec. IIVI we go beyond, investigating in a toy model the full content of the sFKPP equation and its consequences in 
QCD. 

II. UNITARY EVOLUTION OF AMPLITUDES IN THE PARTON MODEL 

In this first section, we recall the physical picture of high energy scattering in the parton model fSec. Ill Ajl insisting 
on the manifest link to reaction-diffusion models. In a second step (Sec. lIIB"|l . following Ref. |24|, we derive the general 
evolution equations in QCD. 

A. High energy scattering as a reaction-diffusion process 

Let us consider the scattering of a hadronic probe off a given target, in the rest frame of the probe and at a fixed 
impact parameter. In the parton model, the target interacts through one of its quantum fluctuations (see Fig. 0. 
The probe effectively "counts" the partons in the current Fock state of the target whose transverse momenta k match 
the one of the probe: the amplitude T(k) for the scattering off this particular partonic configuration is proportional 
to the number of partons n(k). 

In QCD, the wave function of a hadronic object is built up from successive splittings of partons starting from the 
valence structure. As one increases the rapidity Y by boosting the target, the opening of the phase space for parton 
splittings makes the probability for high occupation numbers larger. In the initial stages of the evolution, the parton 
density grows diffusively from these splittings. 

On the other hand, the number of partons in each cell of transverse phase space is effectively limited to a maximum 
number N that depends on the strength of the interaction of the probe with the partons, that is on the relationship 
between T and n. This property is necessary from unitarity, which imposes an upper bound on the amplitude (e.g. 
T(k) < 1 with standard normalizations) which, in turn, results in an upper bound on n(k). This is parton saturation 
and is realized by a recombination process. 



3 



Viewed in this way, scattering in the parton model is a reaction-diffusion process. The rapidity evolution of the 
Fock states of the target hadron is like the time evolution of a set of particles that diffuse in space, multiply and 
recombine so that there is no more than N particles on each site, on the average. We parametrize space by the 
variable x: it will be related to k in QCD. The QCD amplitude T is like the fractional occupation number u = n/N 
in the reaction-diffusion model. 

In the continuum limit, it is known on general grounds |3jj that, up to a change of variable, u obeys the evolution 
equation 

dtu(x,t) — Dd 2 u(x,t) + Xu(x,t) — /iu 2 (x,t) H — = \/2u(x, t) 77 (x, t) , (1) 

v N 

where 77 is a Gaussian white noise, that is, a random function satisfying 

{ri(x, t)) = and (r)(x, t) V (y, t')) = S(t - t')S(x - y) . (2) 

Eq. JD is known as the "Reggeon field theory" equation. It is to be interpreted in the standard Ito way, as the con- 
tinuum limit of a discretized equation. It belongs to the same universality class as the stochasti c Fisher-K olmogorov- 
Petrovsky-Piscounov (sFKPP) equation (which may be obtained by the replacement — > y/u(l — u) in Eq. QJ). 
D and A are parameters that characterize the diffusive growth of the particles. The structure of Eq. (JTJ is very 
transparent: the first two terms represent the growing diffusion of the particles, characterized by the parameters D 
and A, the third term implements the recombination process and tames this growth when maximum occupancy is 
approached. The last term results from the stochastic dynamics, whose effect reflects the finiteness of the number of 
particles. Note that this term gives rise to fluctuations in u{x,t) and, consequently, in the particle number n, that 
are proportional to 5n cx *Jn. This is typical of fluctuations of independent random numbers, and reflects well the 
statistical origin of this term. 

We do not mean that Eq. Q fully represents parton evolution. Indeed, in the parton model, we anticipate that the 
rapidity evolution of parton densities be given, in the linear regime, by the BFKL equation, which is more complicated 
than a second order differential operator. But what we mean is that Eq. describes exactly the critical behavior of 
the QCD amplitude T (~ u), that is its large rapidity (time) and large maximum occupancy number N limit, up to 
the replacement of the relevant parameters. 

All higher order operators beyond the ones appearing in Eq. Q are irrelevant in these limits. As for the linear 
part, this is clear: the large time behavior stems from a saddle point that selects the quadratic approximation to the 
growing diffusion kernel. As a matter of fact, any equation of the form 

d t u(x, t) = K ® u(x, t) — fiu 2 (x, t) H — -=y/2u(x,t) T](x, t) , (3) 

v N 



where K is an appropriate differential or integral kernel encoding the growing diffusion, is believed to belong to that 
universality class and hence to exhibit the same asymptotic solutions. By "appropriate" we mean in practice that the 
phase velocity 

V M = ^-M with if (7) = e^ x (K ® e" 7X ) (4) 
7 

of a wave of wave number 7 propagated by the linear part of Eq. JSJ , must be a function of 7 that has a minimum 
in its domain of analyticity, but we should also stress that a general mathematical definition of the universality class 
of the Reggeon field equation is still not available H^. 1 As for the nonlinear part in Eqs. and ©, terms of the 
form v? for j > 2 would have no influence because the region u ~ 1 in which these terms show up is absorptive [2l| . 
Furthermore, we will see that the time evolution is essentially driven by the region where u is small. Thanks to this 
property, we will not have to bother about the exact way how parton saturation actually occurs in QCD: it is enough 
to know that there is an upper bound on the number of partons by unit of transverse phase space. 

So far, we have focused on the evolution of one (partonic) realization, described by Eq. ©. In the QCD context, 
u ~ T represents the scattering amplitude off one particular Fock state. It is not a physical observable because it is not 
possible to select a particular Fock state of the target in a particle physics experiment. The physical amplitude that 



1 Strictly speaking, the case of the running coupling is not captured by Eq. J3J since K would not be a proper diffusion kernel. However, 
it is possible to treat the running coupling through a generalization of the methods known to treat the sFKPP equation, as was done 
in the mean field case in Ref. |18| . 
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is measured in an experiment is the average over all accessible random (partonic) realizations. Only a few moments of 
T may be measurable: its average, related to total cross sections, and its variance, related to diffractive cross sections, 
namely: 

otot ~ (T) 

a dlff ~ (T 2 ) {T) 2 . W 

The evolution of the moments of T (or of u) may be computed directly from Eq. (or Eq. ©). Taking the average 
of Eq. (J3J, with the help of Eq. (J2J, it is straightforward to obtain the evolution of the latter 

d t (u(x,t)) = K(g> (u(x,t)) - n(u 2 (x,t)} . (6a) 

However, this equation is not closed: it calls for an evolution equation for the correlator (u 2 (x, t)}, that can be derived 
again from Eq. © with the help of Eq. (J2J : 

2 

d t (u(x,t)u(y,t)) = (K x + K y ) ® (u(x,t)u(y,t)) - ^(u 2 (x,t)u{y,t)) - [i{u{x,t)u 2 (y,t)) + —S(x-y)(u(x,t)) , (6b) 

and so on for the higher order correlators, leading eventually to an infinite hierarchy of coupled equations. 

Note that the hierarchical form © was proposed by Parisi and Zhang |3l| to describe a growth process known 
as the Eden model. In the same paper, they also obtained the same hierarchy from the Reggeon action, which had 
already been proved to be related to stochastic processes j33| • Only the boundary terms (last term in the right-hand 
side of Eq. JfjbJ) could not be obtained from the Reggeon action through the Parisi-Zhang procedure: indeed, they 
refer explicitly to the way one counts the objects, as can be seen from the fact that this term is proportional to 2/N, 
and the knowledge of the evolution encoded in the bulk terms of the action they had used is not enough. 

Due to universality of the asymptotics of reaction-diffusion processes, it is straightforward to take over Eq. |J3J| or 
equivalently Eqs. Q6a[) and l|6b|l (and the higher equations in this hierarchy), together with their exact asymptotic 
solutions that will be exhibited in the following sections, to describe high energy QCD. Basically, an elementary 
analysis of the parton splitting process and how the partons interact perturbatively with a probe will be enough. This 
will enable us to identify the variables t and x, the kernel of the linear evolution K and the strength of the noise 1/iV 
that correspond to the QCD parton model. We will turn to this analysis in the next section. 

For Eq. ©, being both explicitly nonlinear and stochastic, it proves a formidable task to find solutions. However, 
there has been much progress recently in this direction |25j |. It was realized that, starting from a localized initial 
condition, the large time solution is a traveling wave moving in the direction of larger x, up to some noise. The 
position X t of that wave may be defined in several ways. For example, one can pick a constant k, and define X t in 
such a way that 

u(X t ,t) = K . (7) 
An alternative definition, useful in particular in the case of discrete models, could be 

oo 



X t = / dxu(x,t) . (8) 
Jo 

In the context of QCD, X t is related to the saturation scale, that is, the point at which parton saturation effects set 
in. 

Outstanding results have been obtained on X t and on the profile of the wave front u(x, t) in the region 1 /N « u < 1. 
An interesting point is that the analytical formulas depend only on some local properties of the characteristic function 
K(j) of the linear diffusion kernel. 



B. The evolution equations of the QCD parton model 

In this section, we review the discussion of Ref. |24j. Our goal is to identify the relevant variables t, x, N and the 
evolution kernel K that correspond to the physical situation of QCD. 

We consider the scattering of a dipole of variable size r (the probe) off a dipole of size r (the target). A natural 
variable that will be used throughout is p=ln(ro/r ). We go to the rest frame of the probe so that the target carries 
all the available rapidity Y . The impact parameter b between the dipoles is fixed. 

At high energy, the quantum fluctuations of the target that interact with the probe are dominated by gluons. It 
proves useful to represent this set of partons by color dipoles |l l| . This is possible in the large- N c limit, where gluons 
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are similar to zero-size qq pairs and non-planar diagrams are suppressed. Subleading-A^ terms do not contribute to 
the evolution of the scattering amplitudes when they are small. However, the dipole approximation breaks down when 
the amplitudes approach their unitarity limits at the same time as nonlinearities set in, but as we will see, that does 
not hamper getting the right asymptotics for the physical quantities that we are going to compute. 

We denote by T(r,Y) the scattering amplitude of the probe off a given partonic realization \u>) of the target, 
obtained after a rapidity evolution Y (the dependence on b is understood). It is a random variable, whose probability 
distribution is related to the distribution of the different Fock state realizations of the target. The values of T(r, Y) 
range between (weak interaction) and 1 (unitarity limit). T(r,Y) will be an essential intermediate quantity in our 
calculations, but as has been explained in Sec. Ill Al it is not an observable. The physical dipole-dipole scattering 
amplitude A(r, Y) is the statistical average over all partonic fluctuations of the target at rapidity Y, i.e. 

A(r,Y)= (T(r,Y)); (9) 

see Eq. (0. 

When T is small, the scattering amplitude off one particular Fock state u>(Y) is the sum of all elementary amplitudes 
with each of the dipoles, i.e. 

T(r,Y)= T el (r,r t ) , (10) 

where i labels the dipoles in the Fock state of the target at the time of the interaction. T e i is the elementary dipole 
interaction and is essentially local in impact parameter. T e i behaves like 

r 2 

Tei(r,n) ~ a 2 s -j- , with r< = min(|r|, \n\) , r> = max(|r|, \n\) (11) 
r > 

when the dipoles overlap, and vanishes otherwise. We have neglected 0(1) factors and logarithms, but these approx- 
imations do not affect the results that we shall obtain, which are largely independent of the details. Eq. (|11|) shows 
that the amplitude T(r, Y) is simply counting the number n(r, Y) of dipoles of size r within a disk of radius r centered 
at the impact parameter of the external dipole: 

T(r,Y)~a 2 s n(r,Y) . (12) 

Note that n{r, Y) can take only discrete values: Indeed, it represents the number of quanta of a particular realization 
of the field of the target. The unitarity bound on T implies that n(r, Y) is also constrained by an upper bound 
N = I /a 2 . Later on, we shall switch to momentum space by the Fourier transformation r <-> k. Of course, the 
relationship 1)12(1 will also hold for the Fourier conjugates. 

We emphasize that the relationship (|12f> is only valid within the framework of the above-mentioned approximations. 
A more accurate treatment of how the counting of the partons is actually realized in QCD would eventually lead to 
more complex equations, see Ref. |27| . We will not enter these complications in this paper since we are only interested 
in the study of the asymptotics that are common with models from statistical physics, for which the counting rule 112|) 
is precise enough. 

The rapidity evolution law for (T(r, Y)) is known and we refer the reader to earlier literature for its derivation 
(e.g. 0). It reads 2 

d Y (T(r,Y)) = [ d 2 z ^ ((T(z,y)) + (T(r-z,Y)) (T(r,Y)) (T(z,Y)T(r- z,Y))) . (13a) 

Eq. H13a|l is not a closed equation for (T): It depends upon the correlator (T(z, Y)T(r—z, Y)). An evolution equation 
for (T(n,Y)T(r 2 , Y)) may be derived in the same way. One gets 

d Y (T(r u Y)T(r 2 ,Y)) = |- J d 2 z -—!j—((T(z,Y)T(r 2 ,Y)) + (T(n- z,Y)T(r 2 ,Y)) - (T(r u Y)T(r 2 ,Y)) 

- (T(z,Y)T(n-z,Y)T(r 2 ,Y))) 

+ T- f d2z K ^ v «T(ri,Y)T(z,Y)) + {T(r u Y)T(r 2 -z,Y)) - {T(n,Y)T(r 2 ,Y)) 
2-k J z A (j 2 -zY 

- {T(n,Y)T(z, Y)T(r 2 -z, F)» , (13b) 



2 The impact parameter dependence could be easily put back in Eq. 1131 . We have omitted it for simplicity and since it is enough for our 
purpose to assume locality of the evolution. 
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which again calls for equations for higher correlators. Eqs. (|13fl turn out to be the first two parts of an infinite 
hierarchy originally derived by Balitsky restricted to dipoles. 

In order to make the correspondence with reaction-diffusion models more concrete at the physical level, it is useful at 
this point to discuss qualitatively the typical shape of T(r, Y) from the evolution. (The reader may find the complete 
derivation in Ref. |24|.) As one takes a new step in rapidity each of the dipoles already present in the wave function 
from the previous step and for which T(ri,Y) 1 may split into two new dipoles. The most probable splittings are 
those in which both child dipoles have a size comparable to the size of the parent dipole: they occur with probability 
one in each interval Ay ~ I/a. Thus the main mechanism for the rise of T(ri,Y) with Y is a growing diffusion 
around the size of the initial dipole r^, that is encoded in the linear part of Eqs. l|13f) . The nonlinear terms appearing 
therein tame this growth in order to satisfy the unitarity limits on (T) . 

In a typical partonic configuration as obtained after a sufficiently large rapidity evolution, the dipoles appear to be 
densely distributed around the size tq (T(r, Y) is large), but they become more rare with decreasing r (or increasing 
p), and for sufficiently large p one meets only rare fluctuations which involve one (or few) dipoles and for which 
T(r) — a 2 s . It is a wave front which with increasing Y progresses towards larger values of p. The position of that front 
is characterized by a momentum scale Q S (Y) called the saturation momentum. It is natural to define it e.g. by the 
value of the inverse dipole size for which the amplitude reaches some predefined number k, i.e. T(r = \/Q S (Y), Y) = n: 
this corresponds to prescription J7J. 

We now proceed to the identification of the QCD evolution equation to the Reggeon field theory equation at the 
technical level. Following Ref. [7J, we perform the transformation 

f°° dr 

T(k,Y)= / -J Q (kr)T(r,Y) (14) 
Jo r 

to momentum space. The main positive outcome of such a transformation is to make the nonlinear terms in Eqs. i|13|) 
local in fc, while obviously the linear part stays an integral kernel. It leads to the formal simplification 

d aY (T(k,Y)) = X (-d L )(T(k,Y)) - (T 2 (k,Y)) 
S a y(T(fc 1> y)T(fc 3 ,y)) = ( X (-d Ll ) + x(~d L2 ))(T(k u Y)T(k 2 ,Y)) 

-{T 2 {k ll Y)T{k 2 ,Y))-{T{k 1 .Y)T 2 (k 2 ,Y)) 



where L = \n(k 2 /A 2 ) and x(l) — 2^(1) — — 7) is the eigenvalue of the BFKL kernel. By x (—Ql), we mean 
the BFKL integro-differential operator that may be defined with the help of the formal series expansion 

X (~d L ) = x(7o)l + x'ho)(-d L - 7ol) + ix"(lo)(-d L - 7o l) 2 + • • • (16) 

for some given 70 between and 1, i.e. for the principal branch of the function y. 

It is clear that, as it stands, the hierarchy H15|) admits the factorized "mean field" solution |3^| 

(T{h, Y) ■ ■ ■ T(k n , y)) = A"- 1 (T(fc 1 , Y)) ■ ■ ■ (T(k n , Y)) (17) 

where A is any constant. However, this is not the physical solution that corresponds to scattering in QCD. The set 
of equations (| 1 51) may be supplemented by boundary terms if one wants to represent explicitly the elementary dipole 
interaction (i.e. the "counting rule"), which is absent in Eq. (compare to Eq. I|6bjl ). By analogy with statistical 
mechanics (see Eq. (|6b(l 'l. this term should be [3ll ] 

2a 2 s 5{Li-L 2 ){T{k u Y)) (18) 

in the second equation. With this term, it is clear that Ijl7|) does no longer solve the hierarchy. Our method of 
derivation does not allow to find directly this boundary term, exactly for the same reasons Parisi and Zhang did 
not get it in their derivation of the hierarchy from the Reggeon action (jHJ . Eq. i|I8|) is nonzero when the transverse 
momenta of the two dipoles in the Fock space of the probe are equal, and thus, may scatter off the same dipole in the 
Fock space of the target: This was neglected in our derivation and has to be reintroduced by hand. 

The evolution equation for T may also be written in the form of a stochastic differential equation, namely, 

d &Y T(k, Y) = X (-d L )T(k, Y) - T 2 (fc, Y) + a s ^2T{k,Y)r](k, Y) (19) 

that gives back the hierarchy l|15|) for the correlators, supplemented with boundary terms such as (|18|l . Such a form 
was proposed for the first time in Refs. |23j | on physical grounds: Indeed, this is the universal asymptotic form of 
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the evolution of any reaction-diffusion type process [30|. Subsequentl y, i t was shown on the technical level how the 
B- JIMWLK formalism can be made consistent with these physics |2g, |27| , and how this term may describe Pomeron 
loops (see also Refs. The source term (TL%1) corresponds exactly to the one found there. 

The identification with the Reggeon field theory equation is now straightforward: 

x — *■ L , t — > aY , 

K -> X (-d L ) , fj,—*l, (20) 
N -> l/o? a . 

Let us discuss the validity of Eq. 1)19(1 when applied to high energy scattering, and the physical expectations. As 
stated before, it is not expected to be exact, but we believe that it captures the essential physics of the QCD parton 
model. Since the linear part drives the front propagation at high rapidity and since it is taken exactly into account 
here (to leading order in a s Y), we expect that the solutions to this equation match the asymptotics of full QCD as far 
as the position of the front is concerned. This means that one should be able to get the asymptotics of the moments 
of the saturation scale Q s for Y large, a s — > 0, and at a fixed impact parameter by solving (|19|l . The latter limitations 
are necessary for the counting rule l(12() to be valid. Similarly, the shape of T in the region 1 /N <C T <C 1 is known 
to be universal and hence should also be obtainable from Eq. ((1911 . 

Of course, strictly speaking, these statements have the status of a conjecture, that will have to be confirmed by 
accurate numerical calculations, and in the long term, by a better mathematical understanding of the solutions of 
equations of the form (|19|l . 

There is still some important difference between the original Reggeon field theory equation Q and the equation 
that we have obtained for high energy QCD. Indeed, T(k, Y) = 1 is not a fixed point of the QCD evolution, unlike 
the case of the original Reggeon field theory equation, so the solutions cannot match in that region. Actually, it is 
known that T(k,Y) ~ \n(Q 2 Jk 2 ) in the region T > 1 35]. But anyway, we would not solve Eq. 1(19(1 around T > 1 
because in the QCD case, the nonlincarity cannot be reduced to a simple quadratic term as in (|19|l . since also color 
structures beyond dipoles are expected to play a role there [29| . That does presumably not influence the asymptotic 
saturation scale, which is completely determined by the linear part of the evolution equation, but it certainly modifies 
the shape of T and of its correlators close to the unitarity limits. 

III. THE MEAN FIELD PICTURE 

The mean field approximation {T'Hk 1 Y)) ~ {T(k, Y)) 2 = A 2 (k, Y) casts Eq. 1(13(1 into a closed form, known as the 
Balitsky-Kovchegov (BK) equation [fjQ: 

d aY A(k, Y) = X {-d L ) A(k, Y) - A 2 (k, Y) . (21) 

A priori, the validity condition for such an approximation is not clear. We will be able to assess it only by studying the 
full problem including fluctuations. It is nevertheless useful to start with such an approximation because it is tractable 
both analytically and numerically, and because the full stochastic solution is nothing else than a perturbation of the 
mean field solution. 



A. General properties of the BK equation 

It has recently been understood that the BK equation belongs to the universality class of the FKPP equation, 
and as such, it admits a family of traveling wave solutions 36]. That is, there exists a function of the rapidity Q S (Y) 
such that 

A(k,Y) = A(L-bxQ 2 s (Y)) (22) 

is a solution to Eq. 121(1 . For this family of solutions, A(£) on the r.h.s. is essentially a slowly varying function of 
£ = L — lnQ^(y) (tending to a constant in the case of the pure FKPP equation) for £ < 0, and is exponentially 
decreasing at large £ 

A® ~ e"T« . (23) 

Q s is the saturation scale introduced before, that is, the transition point between the linear and the saturation regime. 
Note that the general traveling wave property l(22|) is equivalent to geometric scaling, a feature of the data for j*p 
scattering at high energy discovered a few years ago . 
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The crucial point is that the properties of these traveling waves are to a large extent determined by the large-fc tail, 
where the amplitude A is small and thus where Eq. (|21[l may be linearized. This very important property is due to 
the particular propagation mode of the front, which is said to be "pulled" along by its tail. 

The linear part of the BK equation (|21ll has the characteristic function (see Eq. Q) 

«( 7 ) = ^ • (24) 
7 

v (7) is the phase velocity of a wave of wave number 7. It has a minimum at 70 ~ 0.627. 

Starting from a given initial condition A which behaves asymptotically like e - ^, there are two relevant cases. 
Either (3 < 70, in which case the large- K asymptotic solution conserves 7 = /3, or j3 > 70, in which case the wave front 
will converge asymptotically to the shape e -70 ^, that moves at velocity d\nQ 2 s / dY — av(jo). This property can be 
understood in a simple way. The wave packet that corresponds to a physical initial condition, which has at most a 
mild growth as £ — * —00 and that decreases like e - ^ for £ — > +00, contains all waves of wave number 7 ranging from 
-co to (3. At large times, the slowest of these waves will determine the velocity of the wave front. The slowest wave 
is either (3 itself if /3 < 70, or 70 in the opposite case. 

The latter case is the physically relevant one for the BK equation Indeed, color transparency implies that for 
large momenta, i.e. large values of £, the QCD amplitudes behave like e~^. Thus f3 = 1, which is larger than 70. 

The transition from the initial condition at Y = to the asymptotic traveling wave induces corrections to the 
velocity of the front (i.e. to Q s ) and to its shape. The first few orders in a F-expansion of lnQ 2 are completely 
determined by Eq. 



± lnQl{Y) = a^l + Aj^h^ + °( W ■ (25) 

dY 7o 2 7o r 27 2 Va X "( 70 )y3/2 

The first term on the ri ght- hand side was first computed in Ref. Q in the context of the GLR equation. The second 
term was found in Ref. [17J, while the last one was derived for QCD in Ref. Eg. 

The front itself has a rapidity expansion around its asymptotic shape (12. '{1 which reads, to the same order of 
approximation 

/ L2 \ -70 

# ' Y)= iW)j *1>{k,Y), (26) 
where the reduced front ip[k, Y) was computed in Ref. |19| : 



<<i/.-.r. = fv ; " x { -,, iu[/.-W;(r 1] + r-2 + {■>- ^ + l0 ^^f ) - 2 



2 7oX (3) (7o) 1 

3 x"(7o) 3 



i 2 F 2 [l,l;|,3;z 2 ]^z 4 + 6V^(l-iFi [-|, |; z 2 ]) z + 0(l/7F)j . 



(27) 



z = ln(fc 2 IQl{Y))/ \/2ax"{lo)Y is the so-called "leading edge" variable \2^. 

The dominant term (first factor in the right-hand side) is the exponential l|23l) corrected by a linear factor that 
represents absorption, and by a Gaussian whose width increases with rapidity, namely 

Mk , Y) - O, (7. .n (g|y + C,) exp (~ ^™ ) ■ ^ 

Having an explicit Y dependence, this factor clearly violates geometric scaling. One can estimate from Eq. (|28|l how 
many steps in rapidity Ay are needed in order for the exponential shape to set in within a "distance" In k 2 — In Q 2 S 
from the bulk of the front: it corresponds to the point where z approaches 1, i.e. 
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B. Numerical implementation and approximation schemes 

Numerical studies of the BK equation were pioneered in Refs. 0, El El El- These works anticipated qualitatively 
the traveling wave behavior of the large rapidity solution. In this paper, we will focus on the quantitative comparison 
with the most recent analytical results. 

For the purpose of numerical implementation, the following form of the BFKL kernel in Eq. 1211) is well-suited 



X (-d L )T(L) = J + ^dL' 



T(L') - e L - L 'T{L) e L - L 'T(L) 



V4 + e 2 ( L - L ' 



(30) 



There are several known approximations to this kernel. 

In the so-called diffusive approximation, one keeps only terms up to second order in the expansion (|16f) . which 
makes the BFKL kernel local: 

X(7) - x(7o) + x'(7o)(7 " 7o) + ^x"(7o)(7 - 7o) 2 • (31) 

Within this approximation, the BK equation becomes a partial differential equation, which is exactly equivalent to 
the FKPP equation up to a trivial change of variables [lg- Moreover, these terms are enough to yield the exact 
asymptotics of Eqs. Ij25l27|l to that order. For the leading order correction, it is clear because the expansion (pTTf> is 
equivalent to a saddle point approximation. The next-to-leading order only depends on x"(7o)j but this is probably 
an accident, since a priori one would expect x"'(lo) to also play a role. 

One can also take a different point of view, and focus instead on the pole structure of the BFKL kernel, by starting 
from its meromorphic expansion. An accurate representation of the principal branch of the \ function consists in 
keeping only the two poles at 7 — and 7 = 1: 

X ( 7 )~ - + — L_+4(ln2-l) , (32) 
7 1-7 

where the constant is adjusted in such a way that x(i) coincides with its exact value. Within this scheme, x(— 9l) is 
an integral operator that admits the representation 

r+OG c-L 

X {-d L )T{L)~ dL'T(L') + dL' e L '- L T(L') + 4{ln2 - l)T(L) . (33) 

J L J — 00 

We have solved numerically the BK equation with these kernels using several different methods, which has allowed 
cross-checking of the accuracy and reliability. The main algorithm (also used in Ref. (3?|]) is based on Chebyshev 
approximation of the integrand. 3 The integrand is defined on a grid in momentum space given by the extrema of the 
Chebyshev polynomials. The integral is then computed by using a discrete cosine transform, and the ensuing system 
of ordinary differential equations is solved using a fourth-order Runge-Kutta method. This makes the program very 
efficient, and makes it possible to discretize the system on quite a small grid; for most of the results in this paper, a 
grid size of 256 points was sufficient, with the L variable ranging between L m j n = —20 and L max = 138. 

We will now focus on the numerical results and their physical interpretation. 

C. Front formation and wave propagation 

As an initial condition for the evolution, following Ref. [38|, we assume the McLerran-Venugopalan (MV) model 
|39|. defined in coordinate space by 

T(r, Y = 0) = exp(-±r 2 Q 2 s (Y = 0) In [e+l/(r 2 A 2 )] ) , (34) 

with starting scale Q 2 S {Y — 0) set to 1 GeV 2 and A = 200 MeV. The coupling constant a is frozen at 0.2 for all our 
calculations. 



3 The code can be downloaded from http://www.isv.uu.se/~enberg/BK/ 
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Fig. [21 shows the BK evolution of the MV initial condition for rapidities up to Y — 50, displaying a traveling wave 
behavior. The first logarithmic plot clearly shows the steep exponential fall-off (|23[) at very large momenta, while the 
second linear plot shows the mild \a.(Q 2 s / k 2 ) rise towards small momenta. We have checked that the same traveling 
wave front pattern persists for rapidities up to Y — 250. 

To study more precisely the shape of the amplitude, we plot the reduced front tp(k, Y) defined in Eq. (|26|l . which 
enables us to monitor how the asymptotic traveling wave is approached. We display the numerical results of BK 
evolution for different values of the rapidity on Fig. El We compare them to the analytical predictions at both leading 
and next-to-leading order, Eqs. (|28|) and l|27|l . Technically, we start with Eq. I|28|l . and we fit the two constants 
C\ and Ci on the highest rapidity data (Y — 50): we get the values C\ — 0.17 and C2 = 2.9. Then, using this 
determination of the parameters, we compute the reduced front for the lower rapidities, using both the leading order 
formula Eq. I|28|l and the next-to-leading order prediction l|27|l . It turns out that the latter does not differ very much 
from the former. The agreement is not perfect, even when we use the next-to-leading order formula. This reflects the 
fact that rapidities of the order of 50 are still too small for these formulas to apply: the convergence of the asymptotic 
series that gives the shape of the front is only algebraic, and hence very slow. Actually, a better agreement in this 
low rapidity region could be obtained by using our freedom to shift the rapidity Y — > Y + Yq [lij . 

We now turn to the saturation scale. Here, we define it as the value of k for which the amplitude T(k, Y) reaches a 
certain predefined value re (as defined in Eq. (7J|). In Fig. 0] we show an example of amplitudes at several rapidities, 
with the saturation scale marked as a point for each curve for the choice re = 0.01. 

Because the formation of the traveling wave front requires a large rapidity interval, the determination of the 
saturation scale depends on the chosen value of re. We demonstrate this effect in Fig. [31 which shows the analytical 
expressions for the logarithmic derivative of the saturation scale (the velocity of the front) including the first two 
and three terms in the ^-expansion together with the numerical results for three different values of re. Note that the 
smaller the re-value, the longer it takes to reach the asymptotic limit. The larger the re, the closer to the two-term 
analytic result (which needs to be closer to the asymptotic limit to be good) , which at first sight looks quite intriguing 
|40j |. Expressed in another way, different parts of the wave front travel with different velocities, which only become 
equal in the asymptotic Y limit. 

It is clear, however, that the analytical results agrees very well with the numerical results for large Y, where they 
tend to the constant velocity ck s x(7o)/70j corresponding to the exponential increase of the saturation scale. For 
smaller rapidities, we see clearly the effect of subasymptotic effects on the numerical calculation. In particular, the 
dependence of the numerical results on the very definition of the saturation scale is manifest. 

In order to study the agreement of the numerical calculation with Eq. I|25|l at subleading level in more detail we 
define the two quantities 



F(Y) = Y M^(^^-^M + ^) (35 ) 
3 V 2tt \ oY 70 270 Y J 

G{ Y)=Y^ ( dl »S Y) ~ ^ , (36) 



3 V BY 7o 

which, according to the analytical formulas asymptotically behave as F(Y) — > 1 + Q(l/yY) and G(Y) — > — 1 + 
0(1/ VY) (cf. the corresponding expression for the saturation scale in Eq. I|25(l ). In the absence of a third term in 
Eq. I|25|) . F(Y) would instead be close to zero, while G(Y) would be close to — 1. These quantities, extracted from 
the numerical determination of the saturation scale, are plotted in Figs. El and [7\ These plots clearly indicate that 
the third 0(1/ y/Y) -term in the saturation scale is present until rather large rapidities. Furthermore, its influence is 
smaller when the saturation scale is determined "higher up" on the front, where the evolution more quickly approaches 
asymptotics. 

The saturation scale has a non-universal dependence on the initial condition, which is ^-suppressed compared to 
the three terms in the expression (|25|l . To study this dependence we determine the velocity of the evolution from 
three different initial conditions, the MV model, a step function O(L), and a Gaussian centered at k = k . The result, 
see Fig. shows that the three curves converge to the same propagation velocity within a few units in rapidity. The 
corresponding wave fronts are shown in Fig. 

Finally, we would also like to know to what extent the two-pole approximation of the BFKL kernel, represented 
in momentum space by Eq. (|33|l . captures the essentials of the evolution. Fig. 1101 shows a comparison between the 
amplitudes obtained from the two kernels. Evidently they give very similar results. We also compare the velocities in 
Fig. [IT] 

A few comments are in order. There is an obvious difference between the solutions of the FKPP equation, which 
corresponds to the diffusive approximation (|16|) and those of the BK equation. While the FKPP waves are bounded 
by a constant for small k <§; Q s , the BK amplitude diverges like In Q 2 /k 2 , as seen on Fig. [21 This stems from the 
nonlocality of the full BFKL kernel, and can most easily be understood on the simplified form of that kernel in 
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Eq. (|33[1 . Indeed, T(L') = constant does not diagonalize the kernel in the region L — > — oo. Instead, one iteration 
of the kernel generates a factor L s — L = ln(Q 2 /k 2 ) from the first term in Eq. (|33[l . The same phenomenon occurs 
in the large k region, in which the behavior e~ laL gets enhanced by powers of L. However, this does not influence 
the properties of the traveling waves that we are studying here. The reason is that these extra factors are subleading 
with respect to the powers of k 2 that determine the asymptotic behavior of the traveling wave. Our numerical study 
shows an excellent quantitative agreement with the analytical solutions l)25 |l -l|27 |l . 

We note that it is also possible to extract the saturation scale from the numerical results by fitting a functional 
form ln(Q 2 /k 2 ) to the numerical curves for small k |40l |. 

So far, we have focused on a physical initial condition. We considered the MV model, but any initial condition that 
satisfies color transparency at large k 2 would satisfy the condition (3 > 70, and hence, the asymptotic saturation scale 
would obey Eq. (|25Jl and the asymptotic shape of the front at large Y and large k would be T(k) ~ 1/fc 270 . We wish 
however to test further the theoretical framework outlined in Scc lIIIAl bv picking the non-physical initial condition 
T{k) ~ 1/fc, in which case /3 = h < 70- In this case, the asymptotic saturation scale should be 

\nQ 2 s (Y)/dY = a sX {[3)/p ~ 1.1090 • • • (37) 

and the shape of the front should be conserved. This is precisely what is observed in the results of the numerical 
calculation, see Figs. ^] and ED in complete agreement with the theoretical expectations in Sec. IIII Al We note 
however that our results clearly disagree with the claims made in Ref. ^^|, where it was found that the asymptotic 
amplitude goes like T(k) ~ l/k 2l ° at large k in all cases, f3 > 70 and (3 < 70. 



IV. BEYOND THE MEAN FIELD: UNIVERSAL FEATURES FROM A TOY MODEL 

In Sec.[nl we identified the universality class of high energy QCD as that of a stochastic partial differential equation, 
the sFKPP equation. We studied in detail its mean field approximation in Sec. II I II We now aim at investigating the 
effects induced by the stochastic terms neglected in that approximation. 

For this purpose, we will propose a simple toy model that is inspired by the QCD dipole model, and that belongs to 
the same universality class. The advantage of studying a toy model are numerous: First, all universal results can be 
checked and illustrated in a simple and physical way, leaving out the irrelevant complications of QCD. Second, many 
qualitative features, that also hold in the case of QCD, can be investigated. For example, we are going to study in 
detail the relevance of the mean field approximation that was investigated in Sec. IIIII 



A. General results on the sFKPP equation 

Let us start by reviewing the most recent general results that have been obtained regarding universal features of the 
solutions of the sFKPP equation (or equivalently of the Reggeon field theory equation (JTJ). We will then formulate 
our toy model in Sec. II V Bl and illustrate its universal properties outlined here in a numerical solution in Sec. II V CI 

The mean shape of the front in the region where 1/JV < u < 1 was found to be j^H 

u{x,t) sin — e 7lH tJ , (38) 



7770 \ In N 

where 70 minimizes the function 1^(7) of Eq. The large time velocity of the wave front was shown to be 

v = vM - ~^JT • (39) 

This is expected to be the first terms of an asymptotic expansion for large N. Terms of relative order 1/lniV are 
neglected. The scaling of the second moment of the front position is also known from numerical simulations of different 
models [|l| 

o 2 = (X 2 ) (X t ) 2 cx ^ . (40) 

From Eq. (|40|l and (|38|l . we obtain easily the scaling form of the average of u, which corresponds to the physical 
amplitude in the case of QCD |24j : 

(u(x,t))=u\ X ZS^ 1 . (II! 



yJt/]n 3 Nj 
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up to 0(l/y/t) corrections. 

While these results rely mostly on conjectures written down after a numerical study of discrete models |4lj . they 
are believed to have a general validity. Indeed, they are supported both by appealing physical arguments, that we 
will develop below, and by accurate numerical calculations within different specific models. 

In the following, we will propose a toy model that is very similar to QCD and confront the numerical solutions with 
the above-mentioned analytical predictions. 

B. A simple toy model 

We consider a set of particles on a one-dimensional lattice, whose sites are labeled by the coordinate x and separated 
by the distance Ax. These particles may jump from one position to the nearest one on the left or on the right with 
probability pl and p R respectively, and they can multiply with probability A. These elementary processes replace the 
linear diffusive growth of partons as described by the BFKL evolution. In order to enforce a limit on the number of 
particles on the different sites, we impose that each of the n(x, t) particles piled up at site x at time t may die with 
a probability Xn(x, t)/N. This is a way of implementing the equivalent of parton saturation. These rules completely 
define the model. 

This setup is typical of a "reaction-diffusion" process. It could realistically represent, for example, bacterial growth, 
but more generally it represents a wide universality class, to which high energy QCD also belongs, as was explained 
in Sec. II. 

It is straightforward to convert the evolution laws stated above in an equation for the number of particles n{x, t) 
on site x. Between times t and t + At, there are ul{x) particles which jump to the left, n R (x) particles which jump 
to the right, n+(x) which multiply, and n-(x) which disappear. Thus the variation in the particle number reads 

n(x,t+At) - n(x, t) = ~n L (x) - n R (x) - n_(x) + n + {x) + n L (x + Ax) + n R (x- Ax)| t _, t+ At ■ (42) 

^From the rules edicted above, each configuration has a multinomial distribution 

P(|nr,(x), Tl R (x), n+(x). Tl-(x)}) = -. . . . , , . -, \)),' , , , , , . ; rrT 

xp n L -p R -\^(Xn(x,t)/Ny i -(l-p L -p R -X-\n(x,t)/N) An ^^ , (43) 

where An{x,t) = n(x,t) — ni,(x)—n R (x) — n+(x) — n—(x). 

To get a clearer picture of the evolution, one may compute the mean change in the particle number in one time 
step from Eqs. |@2J and (|4*31 : 

(n(x,t + At)\{n(x,t)}) = n(x,t)(l - p L ~pr- \n(x,t)/N) +p L n(x+Ax,t) +p R n(x- Ax, t) + Xn(x,t) . (44) 

Obviously, there are two fixed points of the mean evolution: n(x, t) = and n(x, t) = N. The latter is the maximum 
number of particles allowed on each site, on the average. Let us introduce the fractional occupation number u(x, t) = 
n(x,t)/N. The mean evolution of u(x,t) in one step in time reads 

(u(x,t + At)\{u(x,t)}) = u(x,t) +p L [u(x + Ax 1 t) - u(x,t)] +p R [u(x-Ax,t) - u(x,t)] + Xu(x,t)(l - u(x,t)) . (45) 

Taking the average of both sides of this equation, one sees that it is not a closed equation for (u) . It gets closed only 
after a mean field approximation, obtained through the replacement u — > (it) . 

The phase velocity of a wave of wave number 7 may be obtained from Eq. I|45|) by computing the mean evolution 
of the plane wave e~ ltyX ~ vt ^ in one linear evolution step 

«(7) = ^ln(l + A + PL ( e -^ - 1) + PR (e-< Ax - 1)) . (46) 

This function is the equivalent of 17(7) in Eq. Q in the present context. It is enough to compute the corresponding 
values of 70, u(7o), v "{jo) and to replace them in Eqs. . (|3^|) to get the predictions for the asymptotic shape and 
velocity of the front. 

For the purpose of our numerical study, we choose 

Ax = At = l, p L = PR = 0.1, A = 0.2, (47) 
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QCD 


Toy model 


aY 

«(7) = x{i)h 

L = lnfc 2 /A 2 
L s =lnQ 2 /A 2 
T(k,Y) 

a a 


u(7) = Eq. gnj 

x t 

u(a;, i) 
1/y/N 



TABLE I: Dictionary between QCD and the toy model settings. 

which lead to the following parameters: 

70 = 1.35219 , u(to) = 0.255386 , v"(<y ) = 0.167721 . (48) 

We pick an initial condition that consists of N particles sitting at x — at time t — (u(x = 0, t = 0) = 1). The 
position of the front X t is chosen to be defined by the prescription in Eq. © . 

We evolve that initial condition for different values of N. The result is shown in Fig. El f° r t — 1000. It is a 
front with fluctuations around a mean steady shape. Recall that these fluctuations are statistical fluctuations in the 
particle numbers. We see that their magnitude is consistent with the expectation 

Su ~ y^I/N . (49) 

In view of these results, and in order to be able to go to larger values of N, we choose to replace the stochastic 
evolution (|42|l by its mean field approximation l|44l) in the bins in which the number of particles is larger than say 10 4 , 
together with an appropriate boundary condition on the interface between the stochastic and deterministic calculation. 
According to Eq. (|49|l . this amounts to neglecting effects of order 1% in the concerned bins, which we do not consider 
here and which, empirically, have anyway no sizable influence neither on the form nor on the position of the front. 
Such a hybrid method has been considered before for reaction-diffusion models |42j. 

We may even try to simplify more, keeping the stochasticity only in the foremost bin, along the lines of Ref. |4l| . 
At each new time step, the non-empty bins are updated using a pure mean field approximation. A new bin is filled 
immediately on the right of the rightmost nonempty site, with a number of particles given by a Poisson law of 
parameter 

6 = Nx (u{x,t + At)\{u{x,t)}) . (50) 

where (u(x,t + At)\{u(x,t)}) is given in Eq. I|45|). We will also use this simpler model, to get a feeling on how much 
the solutions depend on the details of the model. For our purpose, it is important to have such indications, since 
eventually, we would like to draw conclusions for QCD from our toy model. 

C. Time dependence of the position of the front 

Starting from our localized initial condition, the first stages of the evolution consist in a diffusive spreading and 
multiplication of the particles, with a gradual filling of the nearby states. Indeed: 

(u(x = 0, At)) = 1 and (u(x = 1, At)) = p R . (51) 

As pr = 0(1) the occupation number of all non-empty sites is very large after this step, and the mean field approxi- 
mation u = (u) is certainly justified. This argument applies also for the few successive steps. Thus, as the mean field 
is justified, we expect that an exponentially decaying wave front be gradually built 

u(x,t) - e^ o{x - Xt) ] (52) 

see Eq. I|23f) . In these early stages, all mean field results can be taken over from Sec. Ill, with the appropriate 
replacements listed in Tab. [I] 

At this point, a comment is in order. The initial condition that we have taken is similar to the McLerran- 
Venugopalan model considered above, which is certainly a good representation of a nucleus or nucleon at very low 
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rapidities. However, in the case of a purely perturbative target, like a single dipole, the equivalent relevant initial 
condition for our toy model would rather be u(x = 0,t = 0) = l/N . Our argument above relies on the large occupation 
numbers and cannot apply here. However, in this case, the initial stages of the evolution consist in building a dense 
state around x = 0, with a steep tail in the direction of x > 0, steeper than the critical front profile e _7oX '. This 
process requires a time of order In JV. Subsequently, the steep profile is converted into the critical one through a mean 
field evolution. 

^From Eq. I|29[) with the appropriate identifications listed in Tab. |U it takes a time of order 

{x-x t y 

for the profile 1">1'I) to diffuse within the distance x — X t from the bulk of the front. During that time, the front velocity 
is given by 

-«■»>>- £ • (54) 

where we have kept only the leading correction in t from Eq. Q25II . The deterministic construction of the front may 
last until the exponential shape has reached the point where u ~ l/N, where it has to stop because of discreteness of 
the number of particles: u(x, t) — 0, l/N, 2/N ■ ■ ■ , and these small discrete numbers can definitely not be interpolated 
by an exponential of the form (|52[1 . /From Eq. i|53|) and from the shape of the asymptotic front Eq. I|52|). we see that 
the latter process stops after the time 

(lnjV/ 7 o) 2 
27 u"(7o) 

at which the velocity of the front is, from Eq. 154[) : 

v(tdiff) = v(jo) 2 . (56) 

cm TV 

If there were no fluctuations in the number of particles, Eq. (|56H would be the asymptotic velocity that takes into 
account the discreteness of u. 

Brunet and Derrida j^J conjectured and checked that, within a specific numerical model, Eq. i|56|) is actually 
the right answer to the stochastic problem. They replaced the full stochastic evolution by a mean field equation 
supplemented by a cutoff that simulates the discreteness in the particle number, and that may be implemented, for 
example, as 

u(x,t + At) = (u(x,t + At)\{u(x,t)})Q[(u(x,t + At)\{u(x,t)}) ~ l/N] , (57) 

where (u(x,t + At)\{u(x,t)}) is given by Eq. (|44|l . They were able to compute the constant c = 6/tt 2 which turns 
Eq. Ij56(l into Eq. I|39|) . and found the profile of the front l|38(l . What Eq. (|56|l actually represents is the mean velocity 
of the front when t 3> In 2 N. The fluctuations that have been neglected so far would lead to a random wandering of 
the position of the front about its mean, but this gives rise to subleading effects in t. Indeed, according to Eq. H4U|) . the 
fluctuations of X t about its mean are proportional to \ft, which is subleading with respect to the mean displacement 
from (jHEJ, proportional instead to t. 

In Figs. 1151 and 1161 we have followed the time evolution of the front velocity of a given realization. The numerical 
calculation is compared to the analytical result for the mean field evolution (|54|l . and the numerical result for the 
modified mean field evolution l|57() . We picked two different values of N for Figs. and ^| respectively. First of all, as 
anticipated, up to a time of order In 2 TV, the stochastic evolution has very few fluctuations, and follows accurately the 
analytical prediction |54jl . and the numerical solution of Eq. (|57J) . The fact that the matching with the leading order 
analytical prediction (|54|) is quantitatively so good is not surprising: indeed, the prescription that we have chosen 
for the measurement of the front velocity, Eq. (|SJ, gives weight to the bulk of the front rather than to the tip, and 
thus would correspond to a value of k of order 1 if the definition J2J were chosen. As demonstrated in our study of 
the BK equation in Sec. Ill, for such values of k, the leading order analytical velocity is an accurate representation 
of the numerics. As soon as we approach t ~ In 2 N , the rise of the front velocity stops. The latter assumes a 
steady mean, consistent with both the analytical calculation 156|) and the numerical solution for the modified mean 
field equation l|57|l . Note however that at this point, the velocity starts to exhibit large short lived fluctuations. In 
addition, sometimes a large fluctuation appears, that lives a time of order In 2 N, and that drives the velocity well 
above its average value expected from mean-field considerations. 



15 



The deterministic building of the front until the time discreteness is felt can be viewed also on a picture of the front 
at different stages of the evolution. This is shown on Fig. El We choose N = 10 10 , and 2 different times: t x = 100, 
which lies within the diffusion time, and ti — 1000, which is much larger. We plot the full front u(x,t) (inset) as well 
as the reduced front e lolyX ~ Xt ^u{x, t). We see that before the diffusion time, the exponential decrease is seen in a very 
limited x-range, whereas once the diffusion time is reached, it extends almost down tou~ l/N. The reduced front for 
t = 100 exhibits the typical Gaussian shape that is expected from a mean field evolution, as in Eq. I|28J) . We have also 
displayed the analytical expectation. The fact that the matching is not perfect in the large x — X t tail is consistent 
with our findings of Sec. Ill in the BK case. Note that the fluctuations are small, and exclusively concentrated in the 
tail of very large x — X t , well ahead of the geometrical scaling zone. For t = 1000, the profile has changed, and looks 
more like an arc of sine as predicted by Eq. ||5SJ|. This time, there are large fluctuations on the edge of the right part 
of the plot. 

Finally, we check quantitatively the validity of the scaling given by Eqs. 1|39|) (or l|5t)|) ) by computing numerically 
v(jo) — (v). The result is displayed on Fig. El Technically, we generate typically 1000 realizations of the evolution of 
the initial condition over a time tf ~ 10 X £<&//, and compute 

{v) = / v(tf)-v(t diff ) \ ^ (5g) 



where the brackets on the r.h.s. stand for an average over the realizations. The result is displayed on Fig. 1181 
together with the analytical expectation. We see that asymptotically, the numerical calculation approaches slowly the 
analytical formula 1)56(1 obtained from the mean field equation with a cutoff (obtained through the Brunet-Derrida 
procedure l|5Tjl and leading to Eqs. Only for very large values of N, the agreement is actually good. For the 

values of N of interest for QCD, typically N = l/a 2 ~ 50 — 100, the analytical calculation is quite far from the 
numerical calculation. However, as can be seen on the figure, the discrepancy amounts essentially to a slowly varying 
function of N. Note that the curve approaches the asymptote from below, which means that the average velocity 
is always larger than the result from the mean field calculation with a cutoff. This can be understood from our 
previous discussion and from Figs. 1151 and 1161 the fluctuations neglected in the Brunet-Derrida procedure pull the 
front ahead of its mean position, and the discrepancy seen on Fig. ITHlmav be traced to the large fluctuations in the 
instantaneous velocity. No more quantitative understanding of these deviations has been achieved up to now. We have 
also plotted on the same figure the result of a calculation in which the evolution equation is mean field everywhere, 
except for the foremost bin, as explained in the introduction to this section. We see an almost perfect matching with 
the fully stochastic model, except for very small values of N: this confirms the observation in j^H that essentially, 
only the stochasticity in the rightmost bin plays a role. However, the discrepancy at small N indicates a breaking of 
universality: the details of the model clearly start to enter. 



D. Variance of the position of the front 



So far, we have followed the evolution of one single realization. Each such realization undergoes a stochastic 
evolution given by Eqs. (|42|l - 143(1 . At each time t, u(x,t) has the universal shape 1(38(1 . up to fluctuations concentrated 
in its tail u(x,t) ~ l/N, see Fig. El As has been checked in Sec. II V CI the average (v) exhibits the i-dependence 
given by Eq. (|56|) . However, since the actual evolution is stochastic, the front position X t undergoes a random walk 
about its mean (X t ), and X t gets a variance er 2 . 

This fact may be understood qualitatively in the following way. Most of the time, the front moves forward in x 
through a deterministic time evolution of its bulk, which has width hiN/jo. But it may happen stochastically that 
several bins get successively filled by particles much ahead of the front. Those particles will then diffuse and multiply, 
and eventually form a new front that will be ahead of the old "deterministic" one. The net effect is a jump in the front 
velocity, one of those that can be seen on Figs. 1151 1161 This phenomenon is of diffusive nature, thus we eventually 
expect the dispersion a 2 of the positions of the front to be proportional to t, as in Eq. (|40p. 

We generate 1000 realizations of the stochastic evolution of our initial condition over the time intervals up to 
t\ = 2000 and £2 = 8000 resp ective ly. The dispersion of the front positions is illustrated in Fig. El We see that 
there is roughly a factor of \rt2Jt1 = 2 in the dispersions between the two considered times, in agreement with 
Eq. (|40|l . Each of these realizations exhibits a universal exponential shape l|52|) in the leading edge region, up to small 
fluctuations concentrated essentially in the tip. However, when taking the average, the curves for different times do 
not superimpose anymore, as seen in the insert of Fig. 1191 Of course, this is due to the dispersion in the position of 
the front gU|). 
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Now we may compute numerically the quantity 

0-2 = / ( x t f - x tdiff - ( x tj - x t d%ss )) 
tf-Uiff \ tf-tdiff 

to check the scaling form 140(1. 1(41(1, The result is displayed in Fig. E01 We see again a good agreement with the 
expectations for large values of N. Note that again, the 1/ln 3 N prediction overestimates the numerical calculation 
for smaller N. We have displayed on the same plot the result of the numerical calculation in the simplified model 
where stochastic effects are concentrated in the rightmost bin. We see again that this model gives the same result 
as the complete one. The discrepancies at large N are statistical fluctuations, due to the fact that the plotted curve 
results from an average over a finite number of realizations (about 1000 here). The higher moments are more and 
more sensitive to such fluctuations. They should disappear when the averaging is done with more realizations. 




E. Consequences for QCD 

So far in this section, we have discussed exclusively the toy model. Let us come back to the case of QCD. The 
analytical results l(38 |) - 1(41 (1 illustrated on the toy model studied in the previous section are readily taken over to QCD 
using Tab. [I] We list them here for completeness. 

The scattering amplitude off a single partonic realization, which is not the physical observable, but is nevertheless 
an important intermediate quantity, reads 

ln(l/ a 2 ) . ( n lQ lnk 2 /Q 2 \ / fc 2 \ ^ 
T(k > Y) ~ Sm I ln(l/o2) ) I Qf J (60) 

up to fluctuations. Note that this formula is only valid in the window 1 <C In k 2 /Q 2 <C In 1 /a 2 . The mean dependence 
of the saturation scale on rapidity reads, asymptotically for large rapidities 

d n n2 \ - X(7o) - 7T 2 7QX"(7o) / ft1 s 
- — (InCJJ = a a n , (61) 

dY X Vs/ 70 21n 2 (l/a 2 ) V ' 

and its variance 




In (l/aj) 



which, together with Eq. 1(60(1 . yields the asymptotic scaling form 



A(k, Y) = (T(k, Y))=A\ hlfc2 (lnQ ' } i ^ ( m ) 



A /aY/ln 3 (l/a 2 ) i 



This last formula quantifies the magnitude of the violations of geometric scaling. Note that Eq. ((61(1 was first obtained 
in Ref. |22j . A violation of geometric scaling was also found there, but the square root was missing in the denominator 
of the argument of A in Eq. 1(63(1 . 

The formulas above are valid as long as aY ^> ln 2 (l/ci! 2 ), and asymptotically for a s <C 1. For uY <C ln 2 (l/a 2 ) 
instead, we have shown that fluctuations do not play a role, and the mean field discussion of Sec. II I II applies. This 
point is particularly clear in Fig. 1171 and Figs. [TBI 1151 The kinematical regime aY -C ln 2 (l/a 2 ) is a window where 
geometric scaling should be at work. 

The results obtained here are exact results of QCD. They are valid for very small values of a s , where perturbation 
theory is justified. So the weak noise limit N 3> 1 or a s <ti 1 that we have considered here is the consistent expansion. 
Actually, as the front of the amplitude travels towards larger values of the momenta, one should take into account the 
running of the coupling, but this is just a (substantial!) technical complication which does not change qualitatively 
the picture. 

Let us discuss now the validity of our results for real world QCD. As can be seen in Figs. 1181 and 1201 the numerical 
calculations get relatively close to the asymptotics at least up to a constant of order 1, for say N > 10 4 , which 
corresponds to a s < 10~ 2 . Realistic attainable values of a s in physics experiments are a s > 10 _1 . This situation 
would rather correspond to a strong noise limit N ~ 1 for which, unfortunately, no predictive theory has yet been 
found. 
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V. CONCLUSION 

In this paper, we focused on the statistical approach to high energy QCD, that we have developed and illustrated. 
We have shown that it provides a very simple picture for high energy scattering, relying directly on the QCD parton 
model. It allows a derivation from first principles of the most sophisticated analytic results that have been obtained 
so far for scattering amplitudes at fixed impact parameter. The simplicity of this approach is the consequence of 
a deep physical connection between the parton model and reaction-diffusion processes. Its main interest is to gain 
a simple physical understanding of the fluctuations present in the extended B-JIMWLK equations, and to enable a 
clear derivation of the universal asymptotics of the scattering amplitudes at high energy. 

We have also solved numerically the QCD evolution equations in the context of a mean field approximation, that 
leads to the BK equation. We have confirmed that the BK equation admits traveling waves, and studied their 
formation starting with different physical initial conditions, with different definitions of the saturation scale. The 
numerical results have been compared to the most recent analytical predictions, and a good agreement has been 
found at all levels, for large enough values of the rapidity. We have also emphasized the role of the different parts of 
the kernel. 

To go beyond the mean field approximation to high energy scattering, we have proposed and solved a toy model 
that captures the essential features of the full QCD evolution. The advantages of studying a toy model are numerous: 
in particular, it is simple enough to allow for light and flexible numerical simulations. Consequently the universal 
properties of the solution can be quite easily studied. We have investigated some of these properties, confirming former 
numerical studies on different models in the same universality class. We have investigated how the asymptotics set 
in. One of the important conclusions of this paper is that in the initial stages of the evolution, the mean field solution 
applies. 

We believe that many features would still deserve to be investigated within such simplified models: for example, 
the effect of the running coupling would be crucial both for phenomenology and for full consistency of the approach, 
that relies on a small-a s asymptotic expansion. This case would correspond to reaction-diffusion in an inhomogeneous 
medium, which in particular allows for a variable maximum of particles per site. 

Obviously, going beyond the leading orders in rapidity and in a s will start to be model dependent (i.e. non- 
universal) at some order, and the complications and specificities of QCD (color, 2-dimensional impact parameter 
space) will show up in the evolution equations. The same is true if one wants to give up the coarse-graining in impact 
parameter |27j. The solutions to the more refined equations, that is the further terms in the asymptotic expansions 
we have considered, will probably be difficult to obtain from the analogy with statistical mechanics that we have 
extensively used in this paper. 



Note added: While this work was being completed, the paper 43] appeared, that also deals with numerical solutions 
of QCD at high energies including fluctuations. The numerical results obtained are complementary to ours since the 
attitude of the author is to focus on values of N realistic in the context of QCD, rather than emphasizing the connection 
to the analytical results obtained on equations in the universality class of the sFKPP equation. 
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Figures 




FIG. 1: The scattering of a probe (here, a virtual photon characterized by the momentum scale k interacting through its qq 
Fock state) off a particular 4-gluon fluctuation of the target dipole. 
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FIG. 2: Evolution of the MV initial condition using the full BK equation I2H . showing T(k, Y) for Y = 0, 5 . . . 50 on a 
logarithmic scale and a linear scale. The Y — curve shows the initial condition. 
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FIG. 3: Fit of the analytical formula for the reduced front to the result from the numerical calculation. Dashed line: leading 
order, Eq. 1271 . Full line: leading+next-to-leading order, Eq. 1281 . Points: result of the numerical integration for various values 
of the rapidity. 
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FIG. 4: Propagating front, with the saturation scale for k — 0.01 marked at each value of rapidity. 
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FIG. 5: Logarithmic derivative of the saturation scale 81n ^^ y ^ obtained from numerical simulations compared to the analytical 
results including two and three terms in the Y-expansion of Eq. 1251 . The saturation scale was obtained by tracking the 
amplitude at heights k = 1.0, 0.01 and 0.0001. The results are shown for two different rapidity ranges, and on the first plot we 
in addition mark the asymptotic velocity a s x(7c)/7c ~ 0.9767. 
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FIG. 6: F(Y) (Eq. £2)) for k = 1.0,0.01 and 0.0001. 
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FIG. 7: G(Y) (Eq. fifty) for k = 1.0,0.01 and 0.0001. 
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FIG. 9: Wave fronts for three different initial conditions: McLerran-Venugopalan, a step function, and a Gaussian (from the 
left to the right). 
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FIG. 10: Comparison between evolution by the lull BK kernel, Eq. 1301 . and the two-pole approximation to the kernel, Eq. 
ITS 
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FIG. 11: Comparison of the logarithmic derivative of the saturation scale 91 "?^ y ' > as obtained from evolution using full BK 
and using the two-pole approximation of the kernel. 
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FIG. 13: Evolution of the shape of the front upon rapidity evolution, starting from an initial condition that behaves like 1/k 
at large k. The rapidity interval ranges from Y — (initial condition) to Y = 100 (from the left to the right). The dotted lines 
represent C'/k. 
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FIG. 14: Numerical integration ol the toy model defined in Eqs. I|42I4^ over 1000 units of time for different values of N. 




FIG. 15: Instantaneous velocity of the front as a function of time (full fluctuating line) from a numerical solution of the stochastic 
evolution equation 11421 for N = 10 6 particles per site. Large dashed curves: analytical mean field solution, asymptotics and 
leading correction. Short dashed curves: numerical solution of the mean field equation supplemented by a cutoff (Eq. JS3) 3 
and analytical asymptotics (Eq. 1391 ') (straight dashed line). 





FIG. 17: Shape of the front after an evolution of the toy model over t\ = 100 (dashed steps) and t2 = 1000 (full steps) steps 
of time. Insert: The front u(x,t) compared to e~ lox . 
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FIG. 18: vo — {v) from the numerical solution of the toy model, see Eq. I|58^ (full curve) compared to the analytical prediction I39|l 
(full straight line). The result from the simplified toy model Eq. 15011 is also displayed (dashed curve). 




FIG. 19: 1000 realizations of the evolution of the toy model between time and t\ = 2000 (left bunch of curves), and = 8000 
(right bunch of curves). Insert: the average front for these two times. 
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FIG. 20: The variance of the front position divided by the time evolution interval, see Eq. I)59|l (full curve). Line: asymptotic 
prediction (1401) (the overall constant being adjusted) . The result from the simplified toy model Eq. 15011 is also displayed (dashed 
curve) . 



